

#include "glodef.h"
#include "errlog.h"
#include "grid3d.h"
#include "fe.h"

#define INTERP_THRESHOLD        LS_ACCURACY * 10

int find_element(const grid3d_t *grid, node_t *nodes, unsigned type,
                 int i, const double *Pt, double **basis);

void transfer_interpolation_vector(const grid3d_t *grid, unsigned type,
                                   node_t *from_nodes, const vector u,
                                   node_t *to_nodes, vector v)
{
    int n, ndof, i, j, Tid, *Tnodes;
    double *pt, *fpt, *basis, val[3];

    assert(grid!=0);
    assert(from_nodes!=0);
    assert(to_nodes!=0);
    assert(u!=0);
    assert(v!=0);
    assert(u!=v);

    if(grid->inverse==NULL)
        error_msg("%s: inverse connectivity is not available.\n", __func__);

    ndof=fe_ndof(grid, type);
    if((get_N_vec(u)!=ndof)||(get_N_vec(v)!=ndof))
        error_msg("%s: wrong dimensions of vectors.\n", __func__);

    if(type!=fe_p1&&type!=fe_bub)
        error_msg("%s: reinterpolation not available for element type %s",
                  __func__, fe_long_names[type]);

    clean_vec(v);

    use_vector(u)
    use_vector(v)
    for(i=0; i<ndof; i++)
    {
        fe_dof_coord(grid, to_nodes, type, i, &pt);
        fe_dof_coord(grid, from_nodes, type, i, &fpt);
        /* check to see if the node moved at all */
        if(dist(pt, fpt)<INTERP_THRESHOLD)
        {
            vX[i]=uX[i];
            vY[i]=uY[i];
            vZ[i]=uZ[i];
        }
        else if(fe_dof_bdry(grid, type, i)>0)
        {
            vX[i]=uX[i];
            vY[i]=uY[i];
            vZ[i]=uZ[i];
        }
        else
        {
            /* find nearest node */
            switch(type)
            {
                case fe_p1:
                    n=i;
                    break;
                case fe_bub:
                    if(i<grid->NN)
                        n=i;
                    else
                        n=grid->elems[i-grid->NN].nodes[0];
                    break;
                case fe_cr:
                default:
                    error_msg("%s: reinterpolation not available for element type %s",
                              __func__, fe_long_names[type]);
            }
            /* find element containing Pt */
            Tid=find_element(grid, from_nodes, type, n, pt, &basis);
            /* interpolate the value of u in from_nodes */
            val[2]=val[1]=val[0]=0.0;
            n=fe_nodes(grid, Tid, type, &Tnodes);
            for(j=0; j<n; j++)
            {
                val[0]+=uX[Tnodes[j]]*basis[j];
                val[1]+=uY[Tnodes[j]]*basis[j];
                val[2]+=uZ[Tnodes[j]]*basis[j];
            }
            /* assign value of v */
            vX[i]=val[0];
            vY[i]=val[1];
            vZ[i]=val[2];
            /* !!! this part is specific to bub elements !!! */
            if(type==fe_bub&&i>=grid->NN)
                for(j=0; j<4; j++)
                {
                    vX[i]-=0.25*vX[Tnodes[j]];
                    vY[i]-=0.25*vY[Tnodes[j]];
                    vZ[i]-=0.25*vZ[Tnodes[j]];
                }
            /* end of part specific to p1 and bub elements */
        }
    }
    done_vector(v)
    done_vector(u)
}

/****************************************************************************
 *
 *
 ****************************************************************************/

int find_element(const grid3d_t *grid, node_t *nodes, unsigned type,
                 int i, const double *Pt, double **basis)
/* this routine attempts to find an element to which the point Pt belongs.
 * i -- index of a node near Pt
 *
 * We check elements containing node i.  After that we check element containing
 * the nodes adjacent to node i.  If we are still unable to locate an element
 * containing P then we trigger an error.
 */
{
    node_inv_t *ninv_i, *ninv;
    int Tid, nn, j;

    if(grid->inverse==NULL)
        error_msg("%s: inverse connectivity is not available.\n", __func__);

    ninv_i=ninv=grid->inverse->ninv+i;
    nn=0;
    while(nn<=ninv_i->NN)
    {
        for(j=0; j<ninv->NEL; j++)
        {
            Tid=ninv->Tids[j];
            if(fe_global_basis(grid, nodes, Tid, type, Pt, basis))
                return Tid;
        }
        ninv=grid->inverse->ninv+ninv_i->nids[nn];
        nn++;
    }
    error_msg("%s: failure.\n", __func__);
}

/****************************************************************************
 *
 *
 ****************************************************************************/
#include "glovar.h"
#include <stdlib.h>
#include <memory.h>

#include "g3dvec.h"

void transfer_taylor_scalar(vector ndisp, unsigned elemP, unsigned elemU)
{
    error_msg("%s is not available because ALE is disabled.\n", __func__);
/*
    fe_matr_p mats[2]={massP, massU};
    // int i;
    scalar tscal, sscal;
    int pN=get_N_scal(p);
    int uN=get_N_vec(u);

    tscal=alloc_scal(pN);
    sscal=alloc_scal(pN);

    // transfer p and phi to the new grid
    // ATTENTION: this only works when p and phi are P1
    if(elemP!=fe_p1)
        error_msg("ALE is only available for pressure in P1.\n");

    size_t foo=grid->NN*sizeof(node_t);
    node_t *nodes=malloc(foo);
    memcpy(nodes, grid->nodes, foo);
    memcpy(grid->nodes, alegrid->old_nodes, foo);

    // node_t *save_nodes = grid->nodes;
    // grid->nodes = alegrid->old_nodes;
    UpdateMeasures(grid);
    fe_matr_assemble(grid, 2, mats);

    PrintLog("Transfering pressure by Taylor expansion.\n");
    test_advection_scalar(grid, elemP, ndisp, p, tscal);
    test_scalar(grid, elemP, p, sscal);
    // for(i=0;i<p->N;i++) tscal->F[i] += sscal->F[i];    
    axpy_scal(tscal, 1.0, sscal);
    fe_matr_solve_scalar(massP, p, tscal, sscal, 0.0);

    //    PrintLog("Transfering phi by Taylor expansion.\n");
    //    test_advection_scalar(grid, elemP, ndisp, phi, tscal);
    //    test_scalar(grid, elemP, phi, sscal);
    //    for(i=0;i<p->N;i++) tscal->F[i] += sscal->F[i];
    //    fe_matr_solve_scalar(massP, phi, tscal, sscal, 0.0);
    
    clean_scal(phi);

    if(!flags.use_ale)
    {
        vector tvec, svec, ndisp_u;
        tvec=alloc_vec(uN);
        svec=alloc_vec(uN);
        ndisp_u=alloc_vec(uN);

        use_vector(ndisp_u)
        use_vector(ndisp)
        memcpy(ndisp_uX, ndispX, grid->NN*sizeof(double));
        memcpy(ndisp_uY, ndispY, grid->NN*sizeof(double));
        memcpy(ndisp_uZ, ndispZ, grid->NN*sizeof(double));
        done_vector(ndisp)
        done_vector(ndisp_u)

        PrintLog("Transfering velocity by Taylor expansion.\n");
        test_advection_vector(grid, elemU, ndisp_u, u, tvec);
        test_vector(grid, elemU, u, svec);
        // for(i=0;i<3*u->N;i++) tvec->X[i] += svec->X[i];
        axpy_vec(tvec, 1.0, svec);
        fe_matr_solve_vector(massU, u, tvec, svec, 0.0);

        free_vec(tvec);
        free_vec(svec);
        free_vec(ndisp_u);
    }
    // grid->nodes = save_nodes;
    memcpy(grid->nodes, nodes, foo);
    free(nodes);
    UpdateMeasures(grid);
    fe_matr_assemble(grid, 2, mats);

    free_scal(tscal);
    free_scal(sscal);
*/
}


